A Counterfactual Simulation Model of Causal Judgment

1 Setup

1.1 Load packages

library("knitr")     # for knitting stuff
library("lme4")      # for linear mixed effects models
library("broom")     # for tidying model results 
library("brms")      # Bayesian regression with Stan
library("corrr")     # for tidy correlation matrix
library("xtable")    # for latex tables
library("jsonlite")  # for reading json files
library("janitor")   # cleans stuff
library("Hmisc")     # bootstrapped confidence intervals
library("tidybayes") # for Bayesian data analysis
library("png")       # adding pngs to images
library("grid")      # functions for dealing with images 
library("egg")       # for geom_custom
library("patchwork") # for figure panels
library("tidyverse") # for wrangling, plotting, etc. 

1.2 Helper functions

# root mean squared error
rmse = function(x, y){
  return(sqrt(mean((x - y)^2)))
}

actual_counterfactual_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = clipindex,
                            y = value,
                            fill = colorindex,
                            group = model)) +
    stat_summary(geom = "bar",
                 fun = "mean",
                 color = "black",
                 position = position_dodge(0.7),
                 width = 0.7) +
    stat_summary(geom = "linerange",
                 fun.data = "mean_cl_boot",
                 position = position_dodge(0.7),
                 size = 1) +
    geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.2) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    facet_grid(index_actual ~ index_counterfactual) +
    geom_text(data = df.text %>%
                drop_na(label),
              mapping = aes(x = x, y = y, label = as.character(label)),
              size = 6,
              position = position_dodge(width = .7),
              hjust = 0.5) +
    coord_cartesian(xlim = c(0.5, 2.5),
                    ylim = c(-5, 105),
                    expand = F,
                    clip = "off") +
    labs(y = ylabel) +
    theme(text = element_text(size = 20),
          legend.position = "none",
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
          panel.spacing.y = unit(0.8, "cm"),
          plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
          axis.title.x = element_blank(),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
  print(p)
}

actual_counterfactual_threeballs_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = ball,
                            y = value,
                            group = model,
                            fill = colorindex)) +
    stat_summary(fun = "mean",
                 geom = "bar",
                 color = "black",
                 position = position_dodge(0.9),
                 width = 0.9) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 color = "black",
                 position = position_dodge(0.9)) +
    geom_point(data = x %>% 
                 mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(jitter.width = 0.3,
                                               jitter.height = 0,
                                               dodge.width = 0.9),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.15) +
    facet_wrap(~clip, ncol = 8) +
    geom_text(data = df.text, 
              mapping = aes(x = x,
                            y = y,
                            label = as.character(label)),
              size = 5, 
              position = position_dodge(width = .9), 
              hjust = 0.5, 
              fontface = "bold") +
    coord_cartesian(xlim = c(0.5, 2.5), 
                    ylim = c(-5, 120), 
                    expand = F,
                    clip = "off") +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    scale_color_manual(values = c("red2", "green2", "black")) +
    labs(y = ylabel) +
    theme(text = element_text(size = 16),
          legend.position = "none",
          axis.title.x = element_blank(),
          panel.spacing.y = unit(c(.3), "cm"),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          strip.text = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
  print(p)
}

# set ggplot theme 
theme_set(theme_classic())

2 Experiment 1: 2 balls

2.1 Clips

include_graphics("../../figures/diagrams/exp1_diagrams.png")

2.2 Read in data

# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv", 
                          header = T,
                          stringsAsFactors = F)

# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# Information about each clip 
df.exp1.clipinfo = tibble(clip = 1:18,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 
                                             1, 1, 0, 1, 1, 1, 1, 1, 1),
                          outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1, 
                                                     0, 1, 1, 0, 0, 1, 0, 1, 1),
                          index_actual = rep(c("actual miss", 
                                               "actual close", 
                                               "actual hit"), each = 6),
                          index_counterfactual = rep(rep(c("counterfactual miss",
                                                           "counterfactual close", 
                                                           "counterfactual hit"),
                                                         each = 2),
                                                     3)) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp1.causal.long = df.exp1.causal$data %>%
  as_tibble() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

Warning: Calling `as_tibble()` on a vector is discouraged, because the behavior is likely to change in the future. Use `tibble::enframe(name = NULL)` instead.
This warning is displayed once per session.

df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
  as_tibble() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

2.2.1 Demographics

# counterfactual condition 
df.exp1.counterfactual %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

   n age_mean age_sd n_female time_mean time_sd
1 41       32    8.8       19       7.5    3.12

# causal condition 
df.exp1.causal %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

   n age_mean age_sd n_female time_mean time_sd
1 41       33   10.3       21      9.65    6.58

2.3 Read in model predictions

2.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]

df.exp1.model = files %>%
  map(~ read.csv(str_c(path, .))) %>%
  reduce(rbind) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp1.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
  group_by(clip) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>%
  ungroup() %>%
  left_join(df.exp1.clipinfo, by = "clip")

# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp1.counterfactual.model = df.exp1.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, ~ sum((.x$prediction*100 - df.exp1.counterfactual.means$rating_mean) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

2.3.2 Causal condition

df.exp1.causal.model = df.exp1.counterfactual.long %>% 
  group_by(clip) %>% 
  summarize(empirical = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.counterfactual.model %>% 
              select(-c(noise, sse)),
            by = "clip") %>% 
  rename(simulation = prediction) %>% 
  mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
         simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))

2.4 Counterfactual condition

2.4.1 Plots

set.seed(1)
df.plot = df.exp1.counterfactual.long %>% 
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.)/2)) %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", "actual close", "actual hit")),
         index_counterfactual = factor(index_counterfactual,
                                       levels = c("counterfactual miss",
                                                  "counterfactual close",
                                                  "counterfactual hit")),
         model = factor(model, levels = c("rating", "prediction"))) %>% 
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12


df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(1:18, each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

actual_counterfactual_plot(x = df.plot, ylabel = "counterfactual judgment")

Warning: Removed 738 rows containing missing values (geom_point).

ggsave("../../figures/plots/exp1_counterfactual_bars.pdf",
       width = 8,
       height = 6)

Warning: Removed 738 rows containing missing values (geom_point).

2.4.2 Stats

2.4.2.1 Model evaluation

df.exp1.counterfactual.means %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(noise = unique(noise), 
            simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 3
  noise simulation_r simulation_rmse
  <dbl>        <dbl>           <dbl>
1     1         0.97            11.0

2.5 Causal condition

2.5.1 Plots

set.seed(1)

model.name = c("empirical")
# model.name = c("simulation")

df.plot = df.exp1.causal.long %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp1.causal.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>%
  pivot_longer(cols = c(rating, simulation, empirical),
               names_to = "model",
               values_to = "value") %>% 
  filter(model %in% c(model.name, "rating")) %>%
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(
    colorindex = ifelse(outcome_actual == 1, 2, 1),
    colorindex = ifelse(model != "rating", 3, colorindex),
    colorindex = as.factor(colorindex),
    index.actual = factor(index_actual, 
                          levels = c("actual miss", "actual close", "actual hit")),
    index.counterfactual = factor(index_counterfactual, 
                                  levels = c("counterfactual miss", 
                                             "counterfactual close", 
                                             "counterfactual hit"))) %>%
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) # swap clips 11 and 12

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(1:18, each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

actual_counterfactual_plot(df.plot, "causal judgment")

Warning: Removed 738 rows containing missing values (geom_point).

ggsave("../../figures/plots/exp1_causal_bars.pdf",
       width = 8,
       height = 6)

Warning: Removed 738 rows containing missing values (geom_point).

2.5.2 Stats

2.5.2.1 Model evaluation

df.exp1.causal.long %>% 
  mutate(rating = abs(rating)) %>% 
  group_by(clip,
           outcome_actual, 
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>% 
  summarize(empirical_r = cor(rating, empirical),
            empirical_rmse = rmse(rating, empirical),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 4
  empirical_r empirical_rmse simulation_r simulation_rmse
        <dbl>          <dbl>        <dbl>           <dbl>
1        0.96           8.57         0.93            15.2

2.5.2.2 Bayesian regression

# random intercepts model 
fit.brm_exp1_causal = brm(formula = rating ~ 1 + index_actual + index_counterfactual + 
                            outcome_actual + (1 + index_actual + index_counterfactual + 
                                                outcome_actual | participant),
                          data = df.exp1.causal.long,
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp1_causal")

2.5.2.3 Specific hypotheses

df.exp1.causal.posterior = df.exp1.causal.long %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp1_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, .value, .draw) %>% 
  pivot_wider(names_from = clip,
              values_from = .value)

func_posterior_difference = function(data, clip1, clip2){
  data %>% 
    mutate(difference = .data[[clip1]] - .data[[clip2]]) %>% 
    pull(difference) %>% 
    mean_hdci() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>% 
    print()  
}

func_posterior_difference(data = df.exp1.causal.posterior,
                          clip1 = "8",
                          clip2 = "13")

            difference
1 (0.26 [-6.51, 7.31])

2.5.3 Tables

fit.brm_exp1_causal %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  xtable() %>%
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:54:37 2020
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
  \toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
intercept & -18.33 & 3.36 & -23.81 & -12.78 \\ 
  index\_actualactualclose & 4.31 & 3.40 & -1.40 & 9.83 \\ 
  index\_actualactualhit & 4.06 & 4.84 & -4.01 & 11.88 \\ 
  index\_counterfactualcounterfactualclose & -17.52 & 3.37 & -22.96 & -12.01 \\ 
  index\_counterfactualcounterfactualhit & -61.37 & 4.38 & -68.54 & -54.30 \\ 
  outcome\_actual & 92.36 & 5.11 & 84.04 & 100.96 \\ 
   \bottomrule
\end{tabular}
\end{table}

3 Experiment 2: Brick & Teleport

3.1 Clips

include_graphics("../../figures/diagrams/exp2_diagrams.png")

3.2 Read in data

df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv",
                                header = T,
                                stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv",
                                        header = T,
                                        stringsAsFactors = F)

# Information about each clip

df.exp2.clipinfo = tibble(clip = 1:20,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 
                                             0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
                          outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 
                                                     1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
                          index_actual = c(rep(c("actual miss", "actual close", "actual hit"),
                                               each = 6),
                                           rep("practice", 2)),
                          index_counterfactual = c(rep(rep(c("counterfactual miss",
                                                             "counterfactual close",
                                                             "counterfactual hit"),
                                                           each = 2), 
                                                       3),
                                                   rep("practice", 2))) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp2.causal_first.long = df.exp2.causal_first$data %>% 
  as_tibble() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("causal", "counterfactual"), each = 20), max(participant)),
         condition = "causal_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>% 
  as_tibble() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("counterfactual", "causal"), each = 20), max(participant)),
         condition = "counterfactual_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
  bind_rows(df.exp2.counterfactual_first.long %>% 
              mutate(participant = participant + 
                       max(df.exp2.causal_first.long$participant)))

3.2.1 Demographics

# counterfactual condition 
df.exp2.counterfactual_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time) %>% 
  bind_rows(df.exp2.causal_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time)) %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

   n age_mean age_sd n_female time_mean time_sd
1 82       35   12.2       45     21.25    5.11

3.3 Read in model predictions

3.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]

df.exp2.model = files %>%
  map(~ read.csv(str_c(path, .))) %>% 
  reduce(rbind) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp2.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp2.counterfactual.means = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  # mutate(rating = abs(rating)) %>% 
  group_by(clip) %>%
  summarize(rating = mean(rating)) %>%
  ungroup() %>%
  left_join(df.exp2.clipinfo, by = "clip")

# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp2.counterfactual.model = df.exp2.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, ~ sum((.x$prediction*100 - df.exp2.counterfactual.means$rating) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

3.3.2 Causal condition

# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
  filter(question == "counterfactual",
    clip <= 18) %>%
  group_by(clip, condition) %>%
  summarize(rating = mean(rating)) %>%
  pivot_wider(names_from = condition,
              values_from = rating) %>% 
  left_join(df.exp2.combined.long %>%
              filter(question == "counterfactual", clip <= 18) %>%
              group_by(clip) %>%
              summarize(combined = mean(rating)),
            by = "clip") %>% 
  left_join(df.exp2.counterfactual.model,
            by = "clip") %>% 
  select(-c(sse, noise)) %>% 
  rename(simulation = prediction) %>% 
  mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined),
            .funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))

3.4 Counterfactual condition

3.4.1 Plots

set.seed(1)

df.plot = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.counterfactual.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")),
         model = factor(model,
                             levels = c("rating", "prediction")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

actual_counterfactual_plot(df.plot, "counterfactual judgment")

# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

3.4.2 Stats

3.4.2.1 Model evaluation

df.exp2.counterfactual.means %>% 
  left_join(df.exp2.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(noise = unique(noise),
            simulation_r = cor(rating, prediction),
            simulation_rmse = rmse(rating, prediction)) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 3
  noise simulation_r simulation_rmse
  <dbl>        <dbl>           <dbl>
1   0.9         0.96            13.5

3.4.2.2 Bayesian regressions

3.4.2.2.1 Check for potential order effects

# random intercepts & slopes model 
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition * 
                                    (index_actual + index_counterfactual) + 
                                    (1 + index_actual + index_counterfactual | participant),
                                  data = df.exp2.combined.long %>% 
                                    na.omit() %>% 
                                    filter(question == "counterfactual"),
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp2_counterfactual")

3.4.3 Tables

fit.brm_exp2_counterfactual %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  xtable() %>%
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:54:41 2020
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
  \toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
intercept & 13.38 & 2.58 & 9.12 & 17.55 \\ 
  conditioncounterfactual\_first & -1.11 & 3.65 & -6.99 & 4.82 \\ 
  index\_actualactualclose & -1.48 & 2.78 & -6.02 & 3.11 \\ 
  index\_actualactualhit & -4.42 & 2.40 & -8.32 & -0.46 \\ 
  index\_counterfactualcounterfactualclose & 37.94 & 3.33 & 32.39 & 43.35 \\ 
  index\_counterfactualcounterfactualhit & 73.80 & 4.17 & 67.00 & 80.57 \\ 
  conditioncounterfactual\_first:index\_actualactualclose & -0.55 & 3.84 & -7.11 & 5.73 \\ 
  conditioncounterfactual\_first:index\_actualactualhit & -0.85 & 3.41 & -6.52 & 4.85 \\ 
  conditioncounterfactual\_first:index\_counterfactualcounterfactualclose & 1.23 & 4.79 & -6.62 & 9.18 \\ 
  conditioncounterfactual\_first:index\_counterfactualcounterfactualhit & -0.02 & 5.92 & -9.86 & 9.55 \\ 
   \bottomrule
\end{tabular}
\end{table}

3.5 Causal condition

3.5.1 Plots

set.seed(1)

# which causal condition? 
# condition.name = c("causal_first")
condition.name = c("counterfactual_first")
# condition.name = c("causal_first", "counterfactual_first")

# which model for the counterfactual simulations 
# model.name = c("simulation")
# model.name = c("causal_first")
# model.name = c("counterfactual_first")
model.name = c("combined")

df.plot = df.exp2.combined.long %>%
  filter(question == "causal", clip <= 18) %>%
  filter(condition %in% condition.name) %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c("rating", model.name),
               names_to = "model",
               values_to = "value") %>% 
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(colorindex = ifelse(outcome_actual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual, 
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(model.name)` instead of `model.name` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

actual_counterfactual_plot(df.plot, "causal judgment")

Warning: Removed 738 rows containing missing values (geom_point).

# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
#        width = 8,
#        height = 6)

3.5.2 Stats

3.5.2.1 Model evaluation

df.exp2.combined.long %>% 
  filter(question == "causal",
         clip <= 18) %>% 
  mutate(rating = abs(rating)) %>% 
  group_by(clip,
           outcome_actual, 
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp2.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>% 
  summarize(empirical_r = cor(rating, combined),
            empirical_rmse = rmse(rating, combined),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 4
  empirical_r empirical_rmse simulation_r simulation_rmse
        <dbl>          <dbl>        <dbl>           <dbl>
1        0.97           8.36         0.94            18.0

3.5.2.2 Bayesian regression

# random intercepts + slopes model 
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition * 
                            (index_actual + index_counterfactual + outcome_actual) +
                            (1 + index_actual + index_counterfactual + 
                               outcome_actual | participant),
                          data = df.exp2.combined.long %>% 
                            filter(question == "causal",
                                   clip <= 18),
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp2_causal")

3.5.2.3 Specific hypotheses

df.exp2.causal.posterior = df.exp2.combined.long %>% 
  filter(question == "causal",
         clip <= 18) %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual, condition) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp2_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, condition, .value, .draw) %>% 
  pivot_wider(names_from = clip,
              values_from = .value)

# difference between clips 13 and 14 versus 8 
df.exp2.causal.posterior %>% 
    mutate(difference = (`13`  + `14`) / 2 - `8`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))

               difference
1 (-0.55 [-12.09, 10.97])

# difference between clips 5 and 6 versus 11 
df.exp2.causal.posterior %>% 
    mutate(difference = (`5`  + `6`) / 2 - `11`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))

             difference
1 (-5.48 [-13.8, 2.73])

3.5.3 Tables

3.5.3.1 Bayesian regression

fit.brm_exp2_causal %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  xtable() %>%
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:54:44 2020
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
  \toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
intercept & -24.07 & 3.50 & -29.80 & -18.27 \\ 
  conditioncounterfactual\_first & 11.47 & 4.91 & 3.29 & 19.40 \\ 
  index\_actualactualclose & 8.17 & 3.44 & 2.39 & 13.62 \\ 
  index\_actualactualhit & 13.42 & 4.90 & 5.42 & 21.30 \\ 
  index\_counterfactualcounterfactualclose & -30.20 & 3.72 & -36.41 & -24.36 \\ 
  index\_counterfactualcounterfactualhit & -50.29 & 4.60 & -57.94 & -42.84 \\ 
  outcome\_actual & 98.53 & 5.40 & 89.53 & 107.57 \\ 
  conditioncounterfactual\_first:index\_actualactualclose & -5.39 & 4.94 & -13.38 & 2.80 \\ 
  conditioncounterfactual\_first:index\_actualactualhit & -16.97 & 7.02 & -28.61 & -5.75 \\ 
  conditioncounterfactual\_first:index\_counterfactualcounterfactualclose & -8.13 & 5.17 & -16.51 & 0.35 \\ 
  conditioncounterfactual\_first:index\_counterfactualcounterfactualhit & -25.51 & 6.44 & -36.53 & -15.18 \\ 
  conditioncounterfactual\_first:outcome\_actual & 10.72 & 7.60 & -1.63 & 23.34 \\ 
   \bottomrule
\end{tabular}
\end{table}

4 Experiment 3: 3 balls

4.1 Clips

include_graphics("../../figures/diagrams/exp3_diagrams.png")

4.2 Read in data

# clipinfo
df.exp3.clipinfo = tibble(clip = 1:32,
                          outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 
                                           0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
                          outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
                          outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
                          outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                                           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                          clipindex = rep(1:2, 16))

# COUNTERFACTUAL JUDGMENTS 
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>% 
  as_tibble() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", "condition"),
                    n()/6),
         participant = rep(1:(n()/6), each = 6)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(condition = factor(condition, levels = 20:21, labels = c("A", "B")),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, condition, gender, age, time, smoothness, difficulty)

# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
  as_tibble() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "rating"), n()/2),
         participant = rep(1:(n()/64), each = 64),
         order = rep(rep(1:32, each = 2), n()/64)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  left_join(df.exp3.counterfactual.demographics %>% 
              select(participant, ball = condition),
            by = "participant") %>% 
  select(participant, ball, clip, everything()) %>% 
  arrange(participant, clip)

# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv",
                          header = T,
                          stringsAsFactors = F)

# demographics 
df.exp3.causal.demographics = df.exp3.causal$extra %>% 
  as_tibble() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", 
                      "counterbalance", "replayType", "experiment"),
                    n()/8),
         participant = rep(1:(n()/8), each = 8)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, gender, age, counterbalance, time, smoothness, difficulty)

# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
  as_tibble() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "x", "y", "z", "rating1", "rating2"), n()/6),
         participant = rep(1:(n()/(32*6)), each = (32*6)),
         order = rep(rep(1:32, each = 6), n()/(32*6))) %>% 
  filter(!name %in% c("x", "y", "z")) %>% 
pivot_wider(names_from = name,
            values_from = value) %>% 
  left_join(df.exp3.causal.demographics %>% 
              select(participant, counterbalance),
            by = "participant") %>% 
  mutate(A = ifelse(counterbalance == 1, rating1, rating2),
         B = ifelse(counterbalance == 1, rating2, rating1)) %>% 
  select(-rating1, -rating2) %>% 
  pivot_longer(cols = c(A, B),
               names_to = "ball",
               values_to = "rating") %>% 
  select(participant, clip, ball, rating, order) %>% 
  arrange(participant, clip, ball)

4.2.1 Demographics

# counterfactual condition 
df.exp3.counterfactual.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

# A tibble: 1 x 6
      n age_mean age_sd n_female time_mean time_sd
  <int>    <dbl>  <dbl>    <int>     <dbl>   <dbl>
1    80       33   10.1       34      18.1    4.63

# causal condition 
df.exp3.causal.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2))

# A tibble: 1 x 6
      n age_mean age_sd n_female time_mean time_sd
  <int>    <dbl>  <dbl>    <int>     <dbl>   <dbl>
1    41       34   10.5       21      21.2    4.96

4.3 Read in model predictions

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]

df.exp3.model = files %>%
  map_dfr(~ read_csv(str_c(path, .))) %>% 
  rename(clip = trial)

4.3.1 Counterfactual condition

# calculate mean counterfactual judgments 
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()

# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp3.model.counterfactual = df.exp3.model %>% 
  group_by(noise) %>% 
  select(clip, contains("whether"), noise) %>% 
  pivot_longer(cols = c(A_whether, B_whether),
               names_to = "ball",
               values_to = "prediction") %>% 
  mutate(ball = str_remove(ball, "_whether")) %>% 
  arrange(noise, clip, ball) %>% 
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, ~ sum((.x$prediction*100 - df.exp3.counterfactual.means$rating_mean) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

4.3.2 Causal condition

# calculate mean causal judgments 
df.exp3.causal.means = df.exp3.causal.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()

df.exp3.model.causal = df.exp3.model %>% 
  # take into account difference-making 
  mutate_at(.vars = vars(A_how:A_robust),
            .funs = list(~ . * A_difference)) %>% 
  mutate_at(.vars = vars(B_how:B_robust),
            .funs = list(~ . * B_difference)) %>% 
  # choose model based on best fit with counterfactual data 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>% 
  pivot_longer(cols = A_difference:B_robust,
               names_to = c("ball", "aspect"),
               names_sep = "_",
               values_to = "value") %>% 
  pivot_wider(names_from = aspect,
              values_from = value) %>% 
  select(clip, ball:robust)

4.3.3 Heuristic model

l.features = fromJSON("data/features.json")

df.features.initial = l.features[["appearance"]] %>%
  cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
  setNames(c(str_c("ball", LETTERS[c(1, 2, 5)], "_appearance"),
             "ballA_velx",
             "ballA_vely",
             "ballB_velx",
             "ballB_vely",
             "ballE_velx",
             "ballE_vely")) %>%
  mutate(clip = 1:nrow(.)) %>%
  pivot_longer(cols = starts_with("ball")) %>%
  separate(name, into = c("ball", "property")) %>%
  pivot_wider(names_from = property, values_from = value) %>% 
  mutate(ball = str_remove_all(ball, pattern = "ball"))

# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")

for (i in 1:32) {
  df.tmp = data.frame()
  for (j in 1:length(ballnames)) {
    ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
    if (ncollisions > 0) {
      tmp = data.frame(
        ball = ballnames[j],
        object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
        time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
        pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
        pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
        post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
        post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
        ncollision = 1:ncollisions
      )
      df.tmp = df.tmp %>%
        rbind(tmp)
    }
  }
  df.features.collisions = df.features.collisions %>%
    rbind(df.tmp %>%
            mutate(clip = i) %>%
            select(clip, ball, everything()))
}

# find end of clip
tmp = df.features.collisions %>%
  group_by(clip) %>%
  filter(ball == "ballE",
         str_detect(object, "goal|Left|Right")) %>%
  group_by(clip) %>%
  mutate(endclip = max(time)) %>%
  select(clip, endclip) %>%
  ungroup() %>%
  distinct()

# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
  left_join(tmp, by = "clip") %>%
  mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
  group_by(clip) %>%
  filter(time <= endclip)

# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
  filter(ball == "E") %>%
  mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
  select(clip, E.moving)

# contact with ball E
tmp.contactE = df.features.collisions %>%
  filter(object == "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  mutate(contactE = 1) %>%
  select(clip, ball, contactE) %>%
  group_by(clip, ball) %>%
  summarize(contactE = any(contactE %>% as.logical()) * 1) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate(contactE = ifelse(is.na(contactE), 0, contactE))

# number of collisions
tmp.ncollisions = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  group_by(clip, ball) %>%
  count() %>%
  rename(ncollision = n)

# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  count(clip, ball, object) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)) %>%
  group_by(clip) %>%
  summarize(AE = any(AE %>% as.logical()) * 1,
            BE = any(BE %>% as.logical()) * 1) %>%
  mutate(A = ifelse(AE == 1 & BE == 0, 1, 0),
         B = ifelse(AE == 0 & BE == 1, 1, 0)) %>%
  select(clip, A, B) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "exclusive")

# impact
func_angle = function(x, y) {
  dot.prod = x %*% y
  norm.x = norm(x, type = "2")
  norm.y = norm(y, type = "2")
  theta = acos(dot.prod / (norm.x * norm.y))
  as.numeric(theta)
}

# impact on ball E 
tmp.impactE = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ball"),
         ball == "ballE") %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff), .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(E.speed.diff = speed.diff,
         E.direction.diff = direction.diff)

# total impact 
tmp.impactTotal = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ballA|ballB")) %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(total.speed.diff = speed.diff,
         total.direction.diff = direction.diff)

# transfer of force 
tmp.transfer = df.features.collisions %>%
  filter(str_detect(ball, "ballA|ballB"),
         str_detect(object, "ball")) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
         AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)) %>%
  mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
  filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
  arrange(clip, time) %>%
  group_by(clip) %>%
  summarize(A = any(!is.na(AE)),
            B = any(!is.na(BE)),
            A = ifelse(any(!is.na(AB)) & any(!is.na(BE)) & max(AB, na.rm = T) < min(BE, na.rm = T),
                       T,
                       A),
            B = ifelse(any(!is.na(AB)) & any(!is.na(AE)) & max(AB, na.rm = T) < min(AE, na.rm = T),
                       T,
                       B)) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "transfer") %>% 
  arrange(clip, ball)

# collect features
df.features = df.features.initial %>%
  filter(ball != "E") %>%
  mutate(ball = factor(ball)) %>%
  mutate(moving = ifelse(velx == 0 & vely == 0, 0, 1),
         speed = abs(velx) + abs(vely)) %>%
  left_join(tmp.contactE) %>%
  left_join(tmp.impactE) %>%
  left_join(tmp.impactTotal) %>%
  left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
  left_join(tmp.movingE) %>%
  left_join(tmp.exclusive) %>%
  select(-c(appearance, velx, vely))

# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])

4.4 Counterfactual condition

4.4.1 Plots

4.4.1.1 Bar graph

set.seed(1)

df.plot = df.exp3.counterfactual.long %>%
  left_join(df.exp3.model.counterfactual %>%
              select(clip, ball, prediction) %>%
              mutate(ball = as.factor(ball)),
            by = c("clip", "ball")) %>%
  left_join(df.exp3.clipinfo,
            by = "clip") %>% 
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = 1,
         colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1,
                             2,
                             colorindex),
         colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1,
                             2,
                             colorindex),
         colorindex = ifelse(model == "prediction", 3, colorindex),
         colorindex = as.factor(colorindex),
         model = factor(model, levels = c("rating", "prediction"))) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>%
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.5,
         colorindex = NA)

actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
#        width = 12,
#        height = 6)

4.4.2 Stats

4.4.2.1 Model evaluation

# best fitting model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.model.counterfactual,
            by = c("clip", "ball")) %>% 
  summarize(simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  mutate_all(~ round(., 2))

Warning: Column `ball` joining factor and character vector, coercing into
character vector

# A tibble: 1 x 2
  simulation_r simulation_rmse
         <dbl>           <dbl>
1         0.86            19.8

# deterministic model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_a, outcome_b) %>% 
              pivot_longer(cols = -clip,
                           names_to = "ball",
                           values_to = "outcome") %>% 
              mutate(ball = factor(ball,
                                   levels = c("outcome_a", "outcome_b"),
                                   labels = c("B", "A"))),
            by = c("clip", "ball")) %>% 
  mutate(outcome = outcome * 100) %>% 
  summarize(simulation_r = cor(rating_mean, outcome),
            simulation_rmse = rmse(rating_mean, outcome)) %>% 
  mutate_all(~ round(., 2))

Warning: Column `ball` joining factors with different levels, coercing to
character vector

# A tibble: 1 x 2
  simulation_r simulation_rmse
         <dbl>           <dbl>
1         0.83            30.3

4.5 Causal condition

4.5.1 Stats

4.5.1.1 Bayesian mixed effects model

df.data = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# using whether-causation from the model 
fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether + (1 + whether | participant),
                                  data = df.data,
                                  save_all_pars = T,
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp3_causal_w")

fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how +
                               (1 + whether + how | participant),
                                   data = df.data,
                                   save_all_pars = T,
                                   cores = 2,
                                   seed = 1,
                                   file = "cache/brm_exp3_causal_wh")

fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient +
                                (1 + whether + how + sufficient | participant),
                                    data = df.data,
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_whs")

# fit.brm_exp3_causal_w %>% summary()
# fit.brm_exp3_causal_wh %>% summary()
# fit.brm_exp3_causal_whs %>% summary()

4.5.1.2 Model comparison

fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
                                            criterion = c("loo", "waic"),
                                            reloo = T)

No problematic observations found. Returning the original 'loo' object.

fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
                                             criterion = c("loo", "waic"),
                                             reloo = T)

No problematic observations found. Returning the original 'loo' object.

fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
                                              criterion = c("loo", "waic"),
                                              reloo = T)

No problematic observations found. Returning the original 'loo' object.

loo_compare(fit.brm_exp3_causal_w,
            fit.brm_exp3_causal_wh,
            fit.brm_exp3_causal_whs)

                        elpd_diff se_diff
fit.brm_exp3_causal_whs    0.0       0.0 
fit.brm_exp3_causal_wh  -125.2      15.6 
fit.brm_exp3_causal_w   -426.2      31.1 

4.5.1.3 Heuristic model

• z-scored predictors and invidiual responses

# Regression based on features
df.regression.features = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  left_join(df.features %>% 
              mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
            by = c("clip", "ball")) %>%
  clean_names() %>% 
  mutate(e_moving = 1 - e_moving)

# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving + 
                                      speed + 
                                      contact_e + 
                                      e_speed_diff + 
                                      e_direction_diff + 
                                      total_speed_diff + 
                                      total_direction_diff + 
                                      transfer + 
                                      e_moving + 
                                      exclusive + (1 | participant),
                                    prior = priors,
                                    data = df.regression.features %>% 
                                      select(-c(clip, ball, order, clipindex,
                                                contains("outcome"))),
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant) 
   Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     8.53      1.22     6.44    11.18 1.00     1163     1936

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept               49.79      1.46    46.91    52.74 1.00     1129
moving                   0.22      0.22     0.01     0.81 1.00     3399
speed                    2.06      0.85     0.44     3.67 1.00     2319
contact_e                0.38      0.35     0.01     1.29 1.00     3201
e_speed_diff             0.12      0.12     0.00     0.46 1.00     4435
e_direction_diff         1.04      0.74     0.04     2.71 1.00     2202
total_speed_diff         2.21      0.95     0.41     4.17 1.00     2208
total_direction_diff     3.99      0.91     2.14     5.68 1.00     2718
transfer                15.59      0.79    14.06    17.14 1.00     3555
e_moving                 0.18      0.17     0.01     0.62 1.00     3296
exclusive                4.39      0.73     2.97     5.81 1.00     3526
                     Tail_ESS
Intercept                1854
moving                   1547
speed                    1266
contact_e                2020
e_speed_diff             2268
e_direction_diff         1685
total_speed_diff         1393
total_direction_diff     1742
transfer                 2992
e_moving                 1738
exclusive                2538

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    33.21      0.46    32.33    34.10 1.00     5313     2845

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

4.5.1.4 Predictions for mean causal judgments

func_fit_data = function(fit){
  fit %>% 
    fitted(newdata = df.exp3.model.causal %>% 
             left_join(df.features %>% 
                         mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
                       by = c("clip", "ball")) %>%
             clean_names() %>% 
             mutate(e_moving = 1 - e_moving),
           re_formula = NA) %>% 
    as_tibble() %>% 
    pull(Estimate)
}

df.exp3.causal.regression = df.exp3.causal.means %>% 
  mutate(w = func_fit_data(fit.brm_exp3_causal_w),
         wh = func_fit_data(fit.brm_exp3_causal_wh),
         whs = func_fit_data(fit.brm_exp3_causal_whs),
         heuristic = func_fit_data(fit.brm_exp3_causal_heuristic))

4.5.1.5 Model evaluation

prediction = "whs"

df.exp3.causal.regression %>% 
  summarize(simulation_r = cor(rating_mean, .[[prediction]]),
            simulation_rmse = rmse(rating_mean, .[[prediction]])) %>% 
  mutate_all(~ round(., 2))

# A tibble: 1 x 2
  simulation_r simulation_rmse
         <dbl>           <dbl>
1         0.87            13.1

4.5.1.6 Individual participant regressions

df.fit = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# priors 
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline = brm(formula = rating ~ 1,
                                              data = df.fit %>% 
                                                filter(participant == 1),
                                              save_all_pars = T,
                                              cores = 2,
                                              seed = 1,
                                              file = str_c("cache/brm_exp3_causal_individual_baseline"))

fit.brm_exp3_causal_individual_w = brm(formula = rating ~ 1 + whether,
                                       data = df.fit %>% 
                                         filter(participant == 1),
                                       prior = priors,
                                       save_all_pars = T,
                                       cores = 2,
                                       seed = 1,
                                       file = str_c("cache/brm_exp3_causal_individual_w"))

fit.brm_exp3_causal_individual_wh = brm(formula = rating ~ 1 + whether + how,
                                        data = df.fit %>% 
                                          filter(participant == 1),
                                        prior = priors,
                                        save_all_pars = T,
                                        cores = 2,
                                        seed = 1,
                                        file = str_c("cache/brm_exp3_causal_individual_wh"))

fit.brm_exp3_causal_individual_whs = brm(formula = rating ~ 1 + whether + how + sufficient,
                                         data = df.fit %>% 
                                           filter(participant == 1),
                                         prior = priors,
                                         save_all_pars = T,
                                         cores = 2,
                                         seed = 1,
                                         file = str_c("cache/brm_exp3_causal_individual_whs"))

# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
  group_by(participant) %>%
  nest() %>%
  ungroup() %>% 
  # fit model for each participant
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.brm_exp3_causal_individual_baseline,
                                          newdata = .x)),
         fit_w = map(.x = data,
                     .f = ~ update(fit.brm_exp3_causal_individual_w,
                                   newdata = .x)),
         fit_wh = map(.x = data,
                      .f = ~ update(fit.brm_exp3_causal_individual_wh,
                                    newdata = .x)),
         fit_whs = map(.x = data,
                       .f = ~ update(fit.brm_exp3_causal_individual_whs,
                                     newdata = .x))) %>% 
  mutate(fit_baseline = map(.x = fit_baseline,
                            ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_w = map(.x = fit_w,
                     ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_wh = map(.x = fit_wh,
                      ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_whs = map(.x = fit_whs,
                       ~ add_criterion(.x, criterion = c("loo", "waic"))),
         r_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ cor(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         r_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ cor(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         r_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ cor(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         rmse_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ rmse(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         rmse_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ rmse(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         rmse_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ rmse(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           w = fit_w,
                                           wh = fit_wh,
                                           whs = fit_whs),
                                 .f = ~ loo_compare(..1, ..2, ..3, ..4)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3", "..4"),
                             labels = c("baseline", "w", "wh", "whs")))

save(list = c("df.exp3.causal.individual_fit"),
     file = "data/exp3_causal_individual_fit.RData")

4.5.1.7 Load individual participant regressions

load("data/exp3_causal_individual_fit.RData")

# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>% 
  count(best_model)

# A tibble: 2 x 2
  best_model     n
  <fct>      <int>
1 wh             2
2 whs           39

4.5.2 Plots

4.5.2.1 Bar graph

set.seed(1)
model_index = "whs"

# model predictions 
model_prediction = list(fit.brm_exp3_causal_w,
                        fit.brm_exp3_causal_wh,
                        fit.brm_exp3_causal_whs) %>%
  map_dfr(~ fitted(., df.exp3.causal.means %>% 
                     left_join(df.exp3.model.causal,
                               by = c("clip", "ball")),
                   re_formula = NA) %>% 
            as_tibble()) %>% 
  mutate(ball = rep(c("A", "B"), n()/2),
         clip = rep(rep(1:32, each = 2), 3),
         model = rep(c("w", "wh", "whs"),
                     each = 64))

df.plot = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_both),
            by = c("clip")) %>% 
  left_join(model_prediction %>% 
              filter(model == model_index) %>% 
              select(model = Estimate, clip, ball),
            by = c("clip", "ball")) %>% 
  pivot_longer(cols = c(rating, model),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = 1,
         colorindex = ifelse(model == "rating" & outcome_both == 0, 1, colorindex),
         colorindex = ifelse(model == "rating" & outcome_both == 1, 2, colorindex),
         colorindex = ifelse(model == "model", 3, colorindex),
         colorindex = as.factor(colorindex),
         model = factor(model, levels = c("rating", "model"))) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>% 
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.5,
         colorindex = NA)

actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
# ggsave("../../figures/plots/exp3_causal_bars.pdf",
#        width = 12,
#        height = 6)

4.5.2.2 Selection

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
  

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(3, 7, 15, 16, 23)) %>%
  group_by(clip, ball) %>%
  summarize(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]) %>%
  left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
            by = c("clip", "ball")) %>%
  ungroup() %>%
  pivot_longer(cols = c(mean, w, wh, whs),
               names_to = "index",
               values_to = "value") %>% 
  mutate_at(.vars = vars(low, high), .funs = ~ ifelse(index == "mean", ., NA)) %>%
  mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
         clip = factor(clip, levels = c(7, 23, 3, 15, 16),
                       labels = c("causal chain", "double prevention", "joint causation",
                                  "overdetermination", "preemption")))

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         value = -10,
         index = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         index = NA,
         value = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = value,
                     group = index,
                     fill = index)) +
  geom_bar(stat = "identity", color = "black",
           position = position_dodge(0.9),
           width = 0.9) +
  geom_errorbar(mapping = aes(ymin = low, ymax = high),
                width = 0,
                size = 1,
                position = position_dodge(0.9)) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -25, ymax = -7) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -25, ymax = -7) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) + 
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_grid(~clip) +
  scale_fill_grey(start = 1,
                  end = 0,
                  labels = c("mean rating  ", expression(CSM[W] ~ " "),
                             expression(CSM[WH] ~ " "),
                             expression(CSM[WHS] ~ " "))) +
  labs(y = "causal responsibility", fill = "") +
  coord_cartesian(clip = "off") + 
  theme_bw() +
  theme(legend.text.align = 0,
        text = element_text(size = 20),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.text = element_text(size = 30),
        strip.text = element_text(size = 30),
        legend.spacing = unit(0.5, "cm"),
        legend.background = element_rect(fill = "transparent"),
        legend.margin = margin(t = 5, unit = "cm"),
        plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))

# ggsave("../../figures/plots/exp3_selection_bars.pdf",
#        width = 20,
#        height = 10,
#        device = cairo_pdf)

4.5.2.3 Scatter plot

func_scatterplot = function(model){
  if(model == "heuristic"){
    xlabel = "Heuristic"
  }else{
    xlabel = bquote(CSM[.(toupper(model))])
  }
  
  l.models = list(w = fit.brm_exp3_causal_w, 
                  wh = fit.brm_exp3_causal_wh,
                  whs = fit.brm_exp3_causal_whs,
                  heuristic = fit.brm_exp3_causal_heuristic)
  
  df.data = df.exp3.causal.means %>% 
    left_join(df.exp3.model.causal, by = c("clip", "ball")) %>% 
    left_join(df.regression.features, by = c("clip", "ball"))
  
  df.plot = l.models[[model]] %>%
    fitted(newdata = df.data,
           re_formula = NA) %>% 
    as_tibble() %>% 
    clean_names() %>%
    bind_cols(df.data)
  
  p = ggplot(data = df.plot,
             mapping = aes(x = estimate, 
                           y = rating_mean)) +
    geom_abline(intercept = 0,
                slope = 1,
                linetype = 2) + 
    geom_smooth(aes(y = estimate,
                    ymin = q2_5,
                    ymax = q97_5),
                stat = "identity",
                color = "black") +
    geom_linerange(size = 0.75, 
                   mapping = aes(ymin = rating_low,
                                 ymax = rating_high),
                   color = "gray50") +
    geom_point(size = 2) + 
    scale_color_manual(values = c("black", "#e41a1c"),
                       guide = F) +
    coord_fixed(xlim = c(0, 100),
                ylim = c(0, 100)) +
    labs(x = xlabel,
         y = "responsibility judgment") +
    annotate(geom = "text",
             label = str_c(
               "r = ", cor(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2) %>% 
                 as.character() %>% 
                 str_sub(start = 2),
               "\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2)),
             x = 5,
             y = 95,
             size = 6,
             hjust = 0) +
    scale_x_continuous(breaks = seq(0, 100, 25)) +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    theme(text = element_text(size = 20))
  return(p)
}

# for creating and saving an individual scatter plot 
# model = "w"
# model = "wh"
# model = "whs"
# model = "heuristic"

# func_scatterplot(model)
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
#        width = 5,
#        height = 5)

list(func_scatterplot(model = "w"),
     func_scatterplot(model = "wh"),
     func_scatterplot(model = "whs"),
     func_scatterplot(model = "heuristic")) %>%
  wrap_plots(ncol = 2)

4.5.2.4 Individual participant variance for a selection of clips (clustered)

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(7, 23, 3, 15, 16)) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
    labels = c("causal chain", "double prevention", "joint causation",
      "overdetermination", "preemption")))

df.cluster = df.exp3.causal.individual_fit %>%
  select(participant, fit_whs) %>%
  mutate(estimates = map(fit_whs, tidy)) %>%
  select(participant, estimates) %>%
  unnest(estimates) %>%
  filter(str_detect(term, "b_"),
         !str_detect(term, "Intercept")) %>%
  mutate(term = str_remove(term, "b_")) %>%
  select(participant, term, estimate) %>%
  pivot_wider(names_from = term,
              values_from = estimate) %>% 
  mutate(group = ifelse(how > whether, "how", "whether"))

df.plot = df.plot %>% 
  left_join(df.cluster %>% 
              select(participant, group),
            by = "participant")

df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         rating = -10,
         group = NA,
         participant = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         group = NA,
         participant = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

# ggplot(df.plot,aes(x=ball,y=rating,group=id,color = best))+
p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = rating,
                     group = participant,
                     shape = group)) +
  geom_line(mapping = aes(linetype = "individual",
                          color = group),
            size = 1,
            alpha = 0.3) +
  geom_point(mapping = aes(color = group),
             size = 1,
             alpha = 0.3) +
  stat_summary(mapping = aes(group = group,
                             color = group,
                             linetype = "mean"),
               fun = "mean",
               geom = "line",
               size = 1.5) +
  stat_summary(data = df.plot %>% filter(group == "how"),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  stat_summary(data = df.plot %>% filter(group == "whether"),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -30, ymax = -10) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -30, ymax = -10) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) +
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_wrap(~clip,
             ncol = 8) +
  coord_cartesian(xlim = c(0.9, 2.1),
                  ylim = c(-5, 105),
                  expand = F,
                  clip = "off") +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     labels = seq(0, 100, 25)) +
  scale_color_brewer(palette = "Set1",
                     guide = "none") +
  labs(y = "causal responsibility rating",
       linetype = "legend",
       shape = "") +
  theme_bw() +
  guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
                                 keywidth = unit(1.2, "cm")),
         shape = guide_legend(override.aes = list(color = c("#e41a1c", "#377eb8"),
                                                  shape = c(19, 19),
                                                  alpha = c(1, 1),
                                                  size = c(5, 5)))) +
  theme(panel.grid = element_blank(),
        text = element_text(size = 20),
        legend.position = c(0.505, 0.23),
        axis.title.x = element_blank(),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 24),
        legend.background = element_rect(fill = "transparent"),
        strip.text = element_text(size = 30),
        axis.text.x = element_blank(),
        legend.key = element_rect(fill = "transparent"),
        legend.box = "horizontal",
        legend.spacing.x = unit(0.1, "cm"),
        panel.spacing.x = unit(1, "cm"),
        plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))

# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
#        plot = p,
#        width = 20,
#        height = 8.5,
#        device = cairo_pdf)

4.5.2.5 Ternary plots for individual regression results

• doesn’t work at the moment: ggtern is incompatible with the most recent ggplot2 update

library("ggtern")
df.plot = df.exp3.causal.individual_fit %>% 
  select(participant, fit_whs) %>% 
  mutate(estimates = map(fit_whs, tidy)) %>% 
  select(participant, estimates) %>% 
  unnest(estimates) %>% 
  filter(str_detect(term, "b_"),
         !str_detect(term, "Intercept")) %>% 
  mutate(term = str_remove(term, "b_")) %>% 
  select(participant, term, estimate) %>% 
  pivot_wider(names_from = term,
              values_from = estimate) %>% 
  mutate_at(vars(how:whether), .funs = list(norm = ~ . / (how + whether + sufficient))) %>% 
  mutate(color = 0,
         color = ifelse(how_norm == max(how_norm), 1, color),
         color = ifelse(whether_norm == max(whether_norm), 2, color),
         color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
         color = factor(color))

ggplot(data = df.plot,
       mapping = aes(x = how,
                     y = sufficient,
                     z = whether,
                     color = color)) + 
  geom_point(alpha = 0.7,
             size = 2,
             show.legend = F) + 
  scale_color_manual(values = c("black", "red", "green", "blue")) + 
  coord_tern() +
  # theme_bw() + 
  theme_showarrows() + 
  theme(text = element_text(size = 20),
                tern.axis.title.T = element_blank(),
                tern.axis.title.L = element_blank(),
                tern.axis.title.R = element_blank())

# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
#        width = 5,
#        height = 5)

4.5.3 Tables

4.5.3.1 Relationship between predictors

df.exp3.model %>% 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
  pivot_longer(cols = A_difference:B_robust,
               names_to = c("ball", "aspect"),
               names_sep = "_",
               values_to = "value") %>% 
  pivot_wider(names_from = aspect,
              values_from = value) %>% 
  select(difference, whether, how, sufficient, robust) %>% 
  correlate() %>% 
  shave() %>% 
  fashion() %>% 
  xtable() %>% 
  print(include.rownames = F,
        booktabs = T)

Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:55:31 2020
\begin{table}[ht]
\centering
\begin{tabular}{llllll}
  \toprule
rowname & difference & whether & how & sufficient & robust \\ 
  \midrule
difference &  &  &  &  &  \\ 
  whether &  .50 &  &  &  &  \\ 
  how &  .79 &  .27 &  &  &  \\ 
  sufficient &  .21 &  .10 &  .36 &  &  \\ 
  robust &  .43 &  .93 &  .24 &  .20 &  \\ 
   \bottomrule
\end{tabular}
\end{table}

4.5.3.2 Population level predictors in the mixed effects models

fit.brm_exp3_causal_w %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(model = "CSM_w") %>% 
  bind_rows(fit.brm_exp3_causal_wh %>% 
              tidy() %>% 
              filter(str_detect(term, "b_")) %>% 
              mutate(model = "CSM_wh")) %>% 
  bind_rows(fit.brm_exp3_causal_whs %>% 
              tidy() %>% 
              filter(str_detect(term, "b_")) %>% 
              mutate(model = "CSM_whs")) %>% 
  mutate(term = tolower(term),
         term = str_remove_all(term, "b_")) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  select(model, everything()) %>% 
  xtable() %>% 
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:55:32 2020
\begin{table}[ht]
\centering
\begin{tabular}{llrrrr}
  \toprule
model & term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
CSM\_w & intercept & 31.68 & 1.61 & 29.02 & 34.37 \\ 
  CSM\_w & whether & 39.30 & 2.16 & 35.76 & 42.82 \\ 
  CSM\_wh & intercept & 10.29 & 1.56 & 7.80 & 12.83 \\ 
  CSM\_wh & whether & 28.40 & 2.68 & 23.96 & 32.70 \\ 
  CSM\_wh & how & 36.07 & 2.60 & 31.80 & 40.27 \\ 
  CSM\_whs & intercept & 10.43 & 1.59 & 7.82 & 13.01 \\ 
  CSM\_whs & whether & 27.93 & 2.60 & 23.67 & 32.18 \\ 
  CSM\_whs & how & 26.47 & 2.91 & 21.76 & 31.35 \\ 
  CSM\_whs & sufficient & 30.90 & 3.06 & 25.77 & 35.96 \\ 
   \bottomrule
\end{tabular}
\end{table}

4.5.3.3 Individual participants regression results

df.exp3.causal.individual_fit %>% 
  select_if(~ (is.numeric(.) | is.factor(.))) %>% 
  mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>% 
  select(participant, everything(), best_model) %>% 
  xtable() %>% 
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:55:32 2020
\begin{table}[ht]
\centering
\begin{tabular}{rrrrrrrl}
  \toprule
participant & r\_w & r\_wh & r\_whs & rmse\_w & rmse\_wh & rmse\_whs & best\_model \\ 
  \midrule
  1 & 0.29 & 0.38 & 0.48 & 30.00 & 29.03 & 27.74 & whs \\ 
    2 & 0.20 & 0.77 & 0.77 & 39.32 & 27.94 & 27.55 & whs \\ 
    3 & 0.37 & 0.78 & 0.79 & 38.90 & 28.00 & 26.95 & whs \\ 
    4 & 0.41 & 0.55 & 0.63 & 43.62 & 40.60 & 38.31 & whs \\ 
    5 & 0.45 & 0.53 & 0.54 & 39.23 & 37.32 & 36.81 & whs \\ 
    6 & 0.18 & 0.54 & 0.54 & 40.84 & 36.32 & 36.07 & whs \\ 
    7 & 0.49 & 0.61 & 0.64 & 32.96 & 29.93 & 29.15 & whs \\ 
    8 & 0.39 & 0.63 & 0.68 & 38.37 & 33.26 & 31.64 & whs \\ 
    9 & 0.40 & 0.72 & 0.76 & 35.86 & 28.34 & 26.39 & whs \\ 
   10 & 0.28 & 0.52 & 0.53 & 29.13 & 26.28 & 25.85 & whs \\ 
   11 & 0.51 & 0.56 & 0.60 & 33.66 & 32.19 & 31.17 & whs \\ 
   12 & 0.56 & 0.57 & 0.70 & 32.78 & 32.08 & 28.84 & whs \\ 
   13 & 0.56 & 0.71 & 0.74 & 29.83 & 25.40 & 24.21 & whs \\ 
   14 & 0.43 & 0.47 & 0.50 & 33.72 & 32.82 & 32.23 & whs \\ 
   15 & 0.19 & 0.35 & 0.37 & 39.77 & 38.29 & 37.94 & whs \\ 
   16 & 0.62 & 0.67 & 0.76 & 40.49 & 37.98 & 34.39 & whs \\ 
   17 & 0.42 & 0.66 & 0.71 & 28.57 & 23.82 & 22.42 & whs \\ 
   18 & 0.31 & 0.35 & 0.39 & 37.78 & 37.09 & 36.57 & whs \\ 
   19 & 0.37 & 0.60 & 0.62 & 33.41 & 29.23 & 28.62 & whs \\ 
   20 & 0.57 & 0.61 & 0.68 & 36.47 & 35.00 & 32.71 & whs \\ 
   21 & 0.43 & 0.50 & 0.57 & 30.95 & 29.61 & 28.32 & whs \\ 
   22 & 0.51 & 0.67 & 0.71 & 33.72 & 29.36 & 28.00 & whs \\ 
   23 & 0.33 & 0.45 & 0.52 & 38.76 & 37.00 & 35.57 & whs \\ 
   24 & 0.46 & 0.63 & 0.69 & 34.01 & 29.96 & 28.22 & whs \\ 
   25 & 0.61 & 0.66 & 0.70 & 30.60 & 28.91 & 27.55 & whs \\ 
   26 & 0.57 & 0.71 & 0.76 & 40.43 & 35.03 & 32.80 & whs \\ 
   27 & 0.46 & 0.52 & 0.66 & 43.12 & 41.49 & 38.12 & whs \\ 
   28 & 0.30 & 0.45 & 0.46 & 29.33 & 27.62 & 27.38 & whs \\ 
   29 & 0.48 & 0.60 & 0.65 & 38.93 & 36.03 & 34.48 & whs \\ 
   30 & 0.18 & 0.70 & 0.70 & 33.06 & 24.95 & 24.71 & whs \\ 
   31 & 0.48 & 0.59 & 0.69 & 37.54 & 34.72 & 31.76 & whs \\ 
   32 & 0.37 & 0.61 & 0.70 & 40.62 & 35.94 & 32.92 & whs \\ 
   33 & 0.50 & 0.56 & 0.64 & 38.64 & 36.79 & 34.85 & whs \\ 
   34 & 0.24 & 0.51 & 0.52 & 44.50 & 40.86 & 40.19 & whs \\ 
   35 & 0.15 & 0.32 & 0.30 & 34.56 & 33.36 & 33.37 & wh \\ 
   36 & 0.19 & 0.84 & 0.82 & 45.38 & 29.58 & 29.70 & wh \\ 
   37 & 0.55 & 0.63 & 0.67 & 36.04 & 33.36 & 32.11 & whs \\ 
   38 & 0.52 & 0.61 & 0.67 & 36.03 & 33.39 & 31.53 & whs \\ 
   39 & 0.47 & 0.77 & 0.77 & 29.61 & 21.75 & 21.58 & whs \\ 
   40 & 0.46 & 0.78 & 0.79 & 32.41 & 23.60 & 22.69 & whs \\ 
   41 & 0.32 & 0.50 & 0.58 & 32.90 & 30.38 & 28.69 & whs \\ 
   \bottomrule
\end{tabular}
\end{table}

4.5.3.4 CSMwhs predictions for selection of cases

df.table = df.exp3.model.causal %>%
  left_join(df.exp3.causal.regression,
            by = c("clip", "ball")) %>% 
  filter(clip %in% c(7, 23, 16, 3, 15)) %>%
  mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic),
            .funs = ~ round(., 2)) %>%
  select(clip, ball, difference, whether, how, sufficient, robust) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ ifelse(. < 0.5,
                             str_c("xmark (", ., ")"),
                             str_c("cmark (", ., ")"))) %>%
  arrange(clip)

xtable(df.table) %>%
  print(include.rownames = FALSE)

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:55:32 2020
\begin{table}[ht]
\centering
\begin{tabular}{lllllll}
  \hline
clip & ball & difference & whether & how & sufficient & robust \\ 
  \hline
7 & A & cmark (1) & xmark (0.34) & cmark (1) & xmark (0) & xmark (0.25) \\ 
  7 & B & cmark (1) & cmark (1) & cmark (1) & cmark (0.67) & cmark (0.6) \\ 
  23 & A & xmark (0.05) & xmark (0) & xmark (0) & xmark (0) & xmark (0) \\ 
  23 & B & cmark (0.91) & cmark (0.79) & xmark (0) & xmark (0) & cmark (0.72) \\ 
  16 & A & cmark (1) & xmark (0.23) & cmark (1) & cmark (1) & xmark (0.35) \\ 
  16 & B & xmark (0) & xmark (0) & xmark (0) & xmark (0) & xmark (0) \\ 
  3 & A & cmark (1) & cmark (0.88) & cmark (1) & xmark (0.12) & cmark (0.76) \\ 
  3 & B & cmark (1) & cmark (0.89) & cmark (1) & xmark (0.11) & cmark (0.75) \\ 
  15 & A & cmark (1) & xmark (0.01) & cmark (1) & cmark (0.99) & xmark (0.1) \\ 
  15 & B & cmark (1) & xmark (0.01) & cmark (1) & cmark (1) & xmark (0.1) \\ 
   \hline
\end{tabular}
\end{table}

4.5.3.5 CSMwhs predictions for all cases

df.table = df.exp3.causal.regression %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball")) %>% 
  left_join(df.exp3.clipinfo %>% 
              select(-clipindex),
            by = c("clip")) %>% 
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ . * 100) %>% 
  mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
  select(clip, ball, contains("outcome"), difference, whether, how, sufficient, robust,
         w, wh, whs, heuristic, rating = rating_mean)

xtable(df.table, digits = 0) %>%
  print(include.rownames = FALSE)

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:55:32 2020
\begin{table}[ht]
\centering
\begin{tabular}{rlrrrrrrrrrrrrrr}
  \hline
clip & ball & outcome\_both & outcome\_a & outcome\_b & outcome\_none & difference & whether & how & sufficient & robust & w & wh & whs & heuristic & rating \\ 
  \hline
1 & A & 0 & 0 & 0 & 0 & 100 & 40 & 100 & 23 & 36 & 47 & 58 & 55 & 57 & 42 \\ 
  1 & B & 0 & 0 & 0 & 0 & 100 & 15 & 100 & 16 & 9 & 38 & 51 & 46 & 54 & 37 \\ 
  2 & A & 0 & 0 & 0 & 0 & 57 & 12 & 0 & 0 & 10 & 36 & 14 & 14 & 25 & 21 \\ 
  2 & B & 0 & 0 & 0 & 0 & 18 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 24 & 19 \\ 
  3 & A & 1 & 0 & 0 & 0 & 100 & 88 & 100 & 12 & 76 & 66 & 71 & 65 & 72 & 76 \\ 
  3 & B & 1 & 0 & 0 & 0 & 100 & 89 & 100 & 11 & 75 & 67 & 72 & 65 & 72 & 75 \\ 
  4 & A & 1 & 0 & 0 & 0 & 100 & 78 & 100 & 4 & 78 & 62 & 69 & 60 & 58 & 63 \\ 
  4 & B & 1 & 0 & 0 & 0 & 100 & 95 & 100 & 15 & 57 & 69 & 73 & 68 & 54 & 78 \\ 
  5 & A & 0 & 0 & 1 & 0 & 100 & 90 & 100 & 0 & 47 & 67 & 72 & 62 & 47 & 71 \\ 
  5 & B & 0 & 0 & 1 & 0 & 100 & 0 & 100 & 0 & 0 & 32 & 46 & 37 & 68 & 22 \\ 
  6 & A & 0 & 0 & 1 & 0 & 100 & 59 & 100 & 16 & 35 & 55 & 63 & 59 & 53 & 73 \\ 
  6 & B & 0 & 0 & 1 & 0 & 100 & 18 & 100 & 6 & 14 & 39 & 52 & 44 & 53 & 22 \\ 
  7 & A & 1 & 0 & 1 & 0 & 100 & 34 & 100 & 0 & 25 & 45 & 56 & 46 & 70 & 59 \\ 
  7 & B & 1 & 0 & 1 & 0 & 100 & 100 & 100 & 67 & 60 & 71 & 75 & 86 & 64 & 79 \\ 
  8 & A & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 25 & 7 \\ 
  8 & B & 1 & 0 & 1 & 0 & 100 & 100 & 100 & 100 & 100 & 71 & 75 & 96 & 84 & 92 \\ 
  9 & A & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 14 & 8 \\ 
  9 & B & 0 & 1 & 0 & 0 & 100 & 100 & 100 & 0 & 100 & 71 & 75 & 65 & 78 & 90 \\ 
  10 & A & 0 & 1 & 0 & 0 & 77 & 18 & 0 & 0 & 22 & 39 & 15 & 16 & 15 & 23 \\ 
  10 & B & 0 & 1 & 0 & 0 & 98 & 79 & 0 & 0 & 63 & 63 & 33 & 33 & 21 & 55 \\ 
  11 & A & 1 & 1 & 0 & 0 & 100 & 70 & 100 & 77 & 68 & 59 & 66 & 80 & 71 & 93 \\ 
  11 & B & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 16 & 4 \\ 
  12 & A & 1 & 1 & 0 & 0 & 100 & 82 & 100 & 74 & 83 & 64 & 70 & 83 & 57 & 77 \\ 
  12 & B & 1 & 1 & 0 & 0 & 100 & 0 & 100 & 12 & 24 & 32 & 46 & 41 & 53 & 37 \\ 
  13 & A & 0 & 1 & 1 & 0 & 67 & 34 & 0 & 0 & 35 & 45 & 20 & 20 & 14 & 8 \\ 
  13 & B & 0 & 1 & 1 & 0 & 70 & 35 & 0 & 0 & 35 & 46 & 20 & 20 & 21 & 64 \\ 
  14 & A & 0 & 1 & 1 & 0 & 97 & 91 & 0 & 0 & 59 & 67 & 36 & 36 & 21 & 22 \\ 
  14 & B & 0 & 1 & 1 & 0 & 91 & 77 & 0 & 0 & 51 & 62 & 32 & 32 & 20 & 18 \\ 
  15 & A & 1 & 1 & 1 & 0 & 100 & 1 & 100 & 99 & 10 & 32 & 47 & 68 & 71 & 76 \\ 
  15 & B & 1 & 1 & 1 & 0 & 100 & 1 & 100 & 100 & 10 & 32 & 47 & 68 & 71 & 76 \\ 
  16 & A & 1 & 1 & 1 & 0 & 100 & 23 & 100 & 100 & 35 & 41 & 53 & 74 & 80 & 92 \\ 
  16 & B & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 14 & 4 \\ 
  17 & A & 0 & 0 & 0 & 1 & 100 & 19 & 100 & 37 & 18 & 39 & 52 & 54 & 66 & 69 \\ 
  17 & B & 0 & 0 & 0 & 1 & 100 & 0 & 100 & 36 & 17 & 32 & 46 & 48 & 65 & 46 \\ 
  18 & A & 0 & 0 & 0 & 1 & 100 & 11 & 100 & 40 & 17 & 36 & 49 & 52 & 55 & 63 \\ 
  18 & B & 0 & 0 & 0 & 1 & 100 & 7 & 100 & 37 & 9 & 34 & 48 & 50 & 56 & 66 \\ 
  19 & A & 1 & 0 & 0 & 1 & 100 & 74 & 100 & 7 & 65 & 61 & 67 & 60 & 55 & 53 \\ 
  19 & B & 1 & 0 & 0 & 1 & 100 & 72 & 100 & 7 & 65 & 60 & 67 & 59 & 55 & 49 \\ 
  20 & A & 1 & 0 & 0 & 1 & 100 & 92 & 100 & 8 & 72 & 68 & 72 & 65 & 57 & 41 \\ 
  20 & B & 1 & 0 & 0 & 1 & 100 & 88 & 100 & 4 & 53 & 66 & 71 & 63 & 56 & 71 \\ 
  21 & A & 0 & 0 & 1 & 1 & 100 & 47 & 100 & 40 & 45 & 50 & 60 & 63 & 58 & 80 \\ 
  21 & B & 0 & 0 & 1 & 1 & 100 & 9 & 100 & 21 & 10 & 35 & 49 & 46 & 59 & 18 \\ 
  22 & A & 0 & 0 & 1 & 1 & 100 & 100 & 100 & 89 & 83 & 71 & 75 & 92 & 48 & 60 \\ 
  22 & B & 0 & 0 & 1 & 1 & 100 & 8 & 100 & 0 & 15 & 35 & 49 & 39 & 53 & 42 \\ 
  23 & A & 1 & 0 & 1 & 1 & 5 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 15 & 3 \\ 
  23 & B & 1 & 0 & 1 & 1 & 91 & 79 & 0 & 0 & 72 & 63 & 33 & 33 & 22 & 39 \\ 
  24 & A & 1 & 0 & 1 & 1 & 100 & 66 & 100 & 4 & 63 & 58 & 65 & 57 & 57 & 44 \\ 
  24 & B & 1 & 0 & 1 & 1 & 100 & 94 & 100 & 22 & 79 & 69 & 73 & 70 & 54 & 73 \\ 
  25 & A & 0 & 1 & 0 & 1 & 100 & 25 & 100 & 21 & 26 & 42 & 53 & 50 & 69 & 43 \\ 
  25 & B & 0 & 1 & 0 & 1 & 100 & 74 & 100 & 54 & 65 & 61 & 67 & 74 & 56 & 73 \\ 
  26 & A & 0 & 1 & 0 & 1 & 100 & 6 & 100 & 3 & 9 & 34 & 48 & 40 & 60 & 39 \\ 
  26 & B & 0 & 1 & 0 & 1 & 100 & 87 & 100 & 35 & 54 & 66 & 71 & 72 & 46 & 69 \\ 
  27 & A & 1 & 1 & 0 & 1 & 100 & 97 & 100 & 52 & 97 & 70 & 74 & 80 & 67 & 80 \\ 
  27 & B & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 17 & 6 \\ 
  28 & A & 1 & 1 & 0 & 1 & 100 & 90 & 100 & 22 & 80 & 67 & 72 & 69 & 79 & 89 \\ 
  28 & B & 1 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 12 & 5 \\ 
  29 & A & 0 & 1 & 1 & 1 & 100 & 58 & 100 & 24 & 44 & 55 & 63 & 61 & 66 & 47 \\ 
  29 & B & 0 & 1 & 1 & 1 & 100 & 63 & 100 & 24 & 38 & 57 & 64 & 62 & 54 & 67 \\ 
  30 & A & 0 & 1 & 1 & 1 & 100 & 57 & 100 & 29 & 49 & 54 & 63 & 62 & 63 & 58 \\ 
  30 & B & 0 & 1 & 1 & 1 & 100 & 46 & 100 & 24 & 39 & 50 & 60 & 57 & 63 & 56 \\ 
  31 & A & 1 & 1 & 1 & 1 & 100 & 2 & 100 & 4 & 3 & 32 & 47 & 39 & 63 & 44 \\ 
  31 & B & 1 & 1 & 1 & 1 & 100 & 4 & 100 & 4 & 4 & 33 & 47 & 39 & 53 & 46 \\ 
  32 & A & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 32 & 10 & 10 & 16 & 5 \\ 
  32 & B & 1 & 1 & 1 & 1 & 100 & 75 & 100 & 66 & 73 & 61 & 68 & 78 & 65 & 71 \\ 
   \hline
\end{tabular}
\end{table}

4.5.3.6 Heuristic model

fit.brm_exp3_causal_heuristic %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  xtable() %>%
  print(include.rownames = F,
        booktabs = T)

% latex table generated in R 3.6.2 by xtable 1.8-4 package
% Sat Mar 21 17:55:32 2020
\begin{table}[ht]
\centering
\begin{tabular}{lrrrr}
  \toprule
term & estimate & std.error & lower 95\% HDI & upper 95\% HDI \\ 
  \midrule
intercept & 49.79 & 1.46 & 47.39 & 52.23 \\ 
  moving & 0.22 & 0.22 & 0.01 & 0.66 \\ 
  speed & 2.06 & 0.85 & 0.64 & 3.43 \\ 
  contact\_e & 0.38 & 0.35 & 0.02 & 1.07 \\ 
  e\_speed\_diff & 0.12 & 0.12 & 0.01 & 0.37 \\ 
  e\_direction\_diff & 1.04 & 0.74 & 0.08 & 2.43 \\ 
  total\_speed\_diff & 2.21 & 0.95 & 0.67 & 3.80 \\ 
  total\_direction\_diff & 3.99 & 0.91 & 2.45 & 5.45 \\ 
  transfer & 15.59 & 0.79 & 14.30 & 16.90 \\ 
  e\_moving & 0.18 & 0.17 & 0.01 & 0.51 \\ 
  exclusive & 4.39 & 0.73 & 3.16 & 5.56 \\ 
   \bottomrule
\end{tabular}
\end{table}

5 Session info

sessionInfo()

R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] forcats_0.5.0   stringr_1.4.0   dplyr_0.8.5     purrr_0.3.3    
 [5] readr_1.3.1     tidyr_1.0.2     tibble_2.1.3    tidyverse_1.3.0
 [9] patchwork_1.0.0 egg_0.4.5       gridExtra_2.3   png_0.1-7      
[13] tidybayes_2.0.1 Hmisc_4.3-1     ggplot2_3.3.0   Formula_1.2-3  
[17] survival_3.1-11 lattice_0.20-40 janitor_1.2.1   jsonlite_1.6.1 
[21] xtable_1.8-4    corrr_0.4.1     brms_2.12.0     Rcpp_1.0.3     
[25] broom_0.5.5     lme4_1.1-21     Matrix_1.2-18   knitr_1.28     

loaded via a namespace (and not attached):
  [1] readxl_1.3.1         backports_1.1.5      plyr_1.8.6          
  [4] igraph_1.2.4.2       splines_3.6.2        svUnit_0.7-12       
  [7] crosstalk_1.1.0.1    rstantools_2.0.0     inline_0.3.15       
 [10] digest_0.6.25        htmltools_0.4.0      rsconnect_0.8.16    
 [13] fansi_0.4.1          magrittr_1.5         checkmate_2.0.0     
 [16] cluster_2.1.0        modelr_0.1.6         matrixStats_0.56.0  
 [19] xts_0.12-0           prettyunits_1.1.1    jpeg_0.1-8.1        
 [22] colorspace_1.4-1     rvest_0.3.5          haven_2.2.0         
 [25] xfun_0.12            callr_3.4.2          crayon_1.3.4        
 [28] zoo_1.8-7            glue_1.3.2           gtable_0.3.0        
 [31] pkgbuild_1.0.6       rstan_2.19.3         abind_1.4-5         
 [34] scales_1.1.0         mvtnorm_1.1-0        DBI_1.1.0           
 [37] miniUI_0.1.1.1       htmlTable_1.13.3     HDInterval_0.2.0    
 [40] foreign_0.8-76       stats4_3.6.2         StanHeaders_2.21.0-1
 [43] DT_0.12              htmlwidgets_1.5.1    httr_1.4.1          
 [46] threejs_0.3.3        arrayhelpers_1.1-0   RColorBrewer_1.1-2  
 [49] ellipsis_0.3.0       acepack_1.4.1        farver_2.0.3        
 [52] pkgconfig_2.0.3      loo_2.2.0            nnet_7.3-13         
 [55] dbplyr_1.4.2         utf8_1.1.4           labeling_0.3        
 [58] tidyselect_1.0.0     rlang_0.4.5          reshape2_1.4.3      
 [61] later_1.0.0          munsell_0.5.0        cellranger_1.1.0    
 [64] tools_3.6.2          cli_2.0.2            generics_0.0.2      
 [67] ggridges_0.5.2       evaluate_0.14        fastmap_1.0.1       
 [70] yaml_2.2.1           processx_3.4.2       fs_1.3.2            
 [73] nlme_3.1-145         mime_0.9             xml2_1.2.5          
 [76] compiler_3.6.2       bayesplot_1.7.1      shinythemes_1.1.2   
 [79] rstudioapi_0.11      reprex_0.3.0         stringi_1.4.6       
 [82] ps_1.3.2             Brobdingnag_1.2-6    nloptr_1.2.2.1      
 [85] markdown_1.1         shinyjs_1.1          vctrs_0.2.4         
 [88] pillar_1.4.3         lifecycle_0.2.0      bridgesampling_1.0-0
 [91] data.table_1.12.8    httpuv_1.5.2         R6_2.4.1            
 [94] latticeExtra_0.6-29  bookdown_0.18        promises_1.1.0      
 [97] boot_1.3-24          colourpicker_1.0     MASS_7.3-51.5       
[100] gtools_3.8.1         assertthat_0.2.1     withr_2.1.2         
[103] shinystan_2.5.0      parallel_3.6.2       hms_0.5.3           
[106] rpart_4.1-15         coda_0.19-3          minqa_1.2.4         
[109] snakecase_0.11.0     rmarkdown_2.1        shiny_1.4.0.2       
[112] lubridate_1.7.4      base64enc_0.1-3      dygraphs_1.1.1.6